function [] = mcsim_icastep()
% 粗对准蒙特卡洛仿真脚本
clear;
clc;
% 加载参数
load('mcsim_ica.mat');
% 修改配置
cfg.fcnPostProc = @mcpp_icastep;
% 修改输入
in.par.refVal(1:3, 2) = 90/3;                                       % 加表零偏，单位μg
in.par.refVal(4:6, 2) = 120/3;                                      % 加速度计标度因数误差，单位ppm
in.par.refVal(7:12, 2) = 73/3;                                      % 加速度计安装误差，单位μrad
in.par.refVal(13:15, 2) = 10/3;                                     % 加速度计二次项系数误差，单位μg/g^2
in.par.refVal(16:18, 2) = [0.012 0.0072 0.036]'/3;                  % 陀螺零偏，单位°/h
in.par.refVal(19:21, 2) = 60/3;                                     % 陀螺标度因数误差，单位ppm
in.par.refVal(22:27, 2) = 73/3;                                     % 陀螺安装误差，单位μrad
in.par.refVal(28, 2) = 2/3;                                         % 加速度计比力斜坡系数，单位μg/h
in.par.refVal(29, 1) = 0.6;                                         % 加速度计比力随机游走系数，单位μg/sqrt(h)
in.par.refVal(30, 1) = 0.005;                                       % 加速度计比速度随机游走系数，单位μg*sqrt(h)
in.par.refVal(31, 1) = 0.6;                                         % 加速度计量化噪声系数，单位m/h
in.par.refVal(32, 2) = 0.12/3;                                      % 陀螺角速度斜坡系数，单位°/h^2
in.par.refVal(33, 1) = 0.005;                                       % 陀螺角速度随机游走系数，单位°/h^(3/2)
in.par.refVal(34, 1) = 0.001;                                       % 陀螺角度随机游走系数，单位°/sqrt(h)
in.par.refVal(35, 1) = 0.05;                                        % 陀螺量化噪声系数，单位″

in.grpSimNum = 200;
% 开始仿真
mcsim(in, cfg);
end

function [out] = mcpp_icastep(in)
% 对准蒙特卡洛仿真后处理函数

dataDir = [pwd '\data\'];
if nargin < 1
    in.grpSimNum = 200;
end
for j = 1:length(in.grpSimNum) % 多组
    close all;
    if in.grpSimNum(j) > 0
%         preparefig('输出值-各轴北向夹角');
        for i = 1:in.grpSimNum(j)
            load([dataDir int2str(j) '_' int2str(i)  '\对准输出分析结果.mat'], 'out_trajGen', 'phiL_CARV2', 'phiL_ICAS', 'out_CARV2', 'out_ICAS');
            % 姿态
            if i == 1
                attError_CARV2 = NaN(length(phiL_CARV2), 3, in.grpSimNum(j));
                attError_ICAS = NaN(length(phiL_ICAS), 3, in.grpSimNum(j));
            end
            attError_CARV2(:, :, i) = phiL_CARV2./pi*180*3600; % 角秒
            attError_ICAS(:, :, i) = phiL_ICAS./pi*180*3600; % 角秒
%             for k = 1:3
%                 subplot(3, 2, k*2-1);
%                 plot(out_trajGen.t(2:end, 1), out_ICAS.azm(:, k)/pi*180);
%                 hold on;
%             end
%             for k = 1:3
%                 subplot(3, 2, k*2);
%                 plot(out_trajGen.t(2:end, 1), out_CARV2.azm(:, k)/pi*180);
%                 hold on;
%             end
        end
        % 绘制方位误差分布图
        preparefig('解析式各次仿真方位误差离散图');
        stem(reshape(attError_CARV2(end, end, :), 1, in.grpSimNum(j)));
        xlabel('仿真次数 / 次'); ylabel('方位误差 / ″');
        preparefig('解析式方位误差频数直方图');
        hist(reshape(attError_CARV2(end, end, :), 1, in.grpSimNum(j)), 20);
        xlabel('方位误差 / ″'); ylabel('频数 / 次');
        
        preparefig('惯性系各次仿真方位误差离散图');
        stem(reshape(attError_ICAS(end, end, :), 1, in.grpSimNum(j)));
        xlabel('仿真次数 / 次'); ylabel('方位误差 / ″');
        preparefig('惯性系方位误差频数直方图');
        hist(reshape(attError_ICAS(end, end, :), 1, in.grpSimNum(j)), 20);
        xlabel('方位误差 / ″'); ylabel('频数 / 次');
        % 绘制收敛图
        t = out_trajGen.t(2:end, 1);
        preparefig('phi_x_CARV2收敛图');
        for i = 1:in.grpSimNum(j)
            plot(t, attError_CARV2(:, 1, i)); hold on;
        end
        grid on; xlabel('time / s'); ylabel('\phi_x^L / (\prime\prime)');
        preparefig('phi_y_CARV2收敛图');
        for i = 1:in.grpSimNum(j)
            plot(t, attError_CARV2(:, 2, i)); hold on;
        end
        grid on; xlabel('time / s'); ylabel('\phi_y^L / (\prime\prime)');
        preparefig('phi_z_CARV2收敛图');
        for i = 1:in.grpSimNum(j)
            plot(t, attError_CARV2(:, 3, i)); hold on;
        end
        grid on; xlabel('time / s'); ylabel('\phi_z^L / (\prime\prime)');
        
        preparefig('phi_x_ICAS收敛图');
        for i = 1:in.grpSimNum(j)
            plot(t, attError_ICAS(:, 1, i)); hold on;
        end
        grid on; xlabel('time / s'); ylabel('\phi_x^L / (\prime\prime)');
        preparefig('phi_y_ICAS收敛图');
        for i = 1:in.grpSimNum(j)
            plot(t, attError_ICAS(:, 2, i)); hold on;
        end
        grid on; xlabel('time / s'); ylabel('\phi_y^L / (\prime\prime)');
        preparefig('phi_z_ICAS收敛图');
        for i = 1:in.grpSimNum(j)
            plot(t, attError_ICAS(:, 3, i)); hold on;
        end
        grid on; xlabel('time / s'); ylabel('\phi_z^L / (\prime\prime)');
        % RMS统计
        out.phiRMS_CARV2 = 3*rms(attError_CARV2, 3);
        out.phiRMS_ICAS = 3*rms(attError_ICAS, 3);
        % 绘制姿态误差
        t = out_trajGen.t(2:end, 1);
        preparefig('phi_x_CARV2');
        plot(t, out.phiRMS_CARV2(:, 1), 'r-', 'LineWidth', 2);
        grid on; xlabel('time / s'); ylabel('\phi_x^L / (\prime\prime)');
        preparefig('phi_y_CARV2');
        plot(t, out.phiRMS_CARV2(:, 2), 'r-', 'LineWidth', 2);
        grid on; xlabel('time / s'); ylabel('\phi_y^L / (\prime\prime)');
        preparefig('phi_z_CARV2');
        plot(t, out.phiRMS_CARV2(:, 3), 'r-', 'LineWidth', 2);
        grid on; xlabel('time / s'); ylabel('\phi_z^L / (\prime\prime)');
        
        preparefig('phi_x_ICAS');
        plot(t, out.phiRMS_ICAS(:, 1), 'r-', 'LineWidth', 2);
        grid on; xlabel('time / s'); ylabel('\phi_x^L / (\prime\prime)');
        preparefig('phi_y_ICAS');
        plot(t, out.phiRMS_ICAS(:, 2), 'r-', 'LineWidth', 2);
        grid on; xlabel('time / s'); ylabel('\phi_y^L / (\prime\prime)');
        preparefig('phi_z_ICAS');
        plot(t, out.phiRMS_ICAS(:, 3), 'r-', 'LineWidth', 2);
        grid on; xlabel('time / s'); ylabel('\phi_z^L / (\prime\prime)');
        
        %% 生成报告
        saveDataDir = [dataDir '第' int2str(j) '组蒙特卡洛计算结果\'];
        if ~isfolder(saveDataDir)
            mkdir(saveDataDir);
        end
        diary([saveDataDir int2str(in.grpSimNum(j)) '次蒙特卡洛计算结果.txt']);
        diary on;
        % 输出时间
        disp(['时间：' datestr(now)]);
        % 输出最终姿态误差
        disp('解析式对准姿态误差3*RMS（单位 角秒）');
        disp(out.phiRMS_CARV2(end, :));
        disp('惯性系对准姿态误差3*RMS（单位 角秒）');
        disp(out.phiRMS_ICAS(end, :));
        diary off;
        saveallfigs(saveDataDir);
    end
end
end